Hopping dynamics for localized Lyapunov vectors in many-hard-disk systems 
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The dynamics of the localized region of the Lyapunov vector for the largest Lyapunov exponent is 
discussed in quasi-one-dimensional hard-disk systems at low density. We introduce a hopping rate to 
quantitatively describe the movement of the localized region of this Lyapunov vector, and show that 
it is a decreasing function of hopping distance, implying spatial correlation of the localized regions. 
This behavior is explained quantitatively by a brick accumulation model derived from hard-disk 
dynamics in the low density limit, in which hopping of the localized Lyapunov vector is represented 
as the movement of the highest brick position. We also give an analytical expression for the hopping 
rate, which is obtained us a sum of probability distributions for brick height configurations between 
two separated highest brick sites. The results of these simple models are in good agreement with 
the simulation results for hard-disk systems. 

PACS numbers: 05.45.Jn, 05.45.Pq, 02.70.Ns, 05.20.Jj 



I. INTRODUCTION 



The dynamical instability and many-body nature play 
essential roles in the justification of a statistical treat- 
ment for deterministic dynamical systems. The dynam- 
ical instability is described as a rapid expansion of the 
difference between two nearby trajectories, namely the 
Lyapunov vector, and the system is called chaotic if at 
least one exponential rate (Lyapunov exponent) of diver- 
gence or contraction of the amplitude of Lyapunov vector 
is positive. The Lyapunov exponent A is defined for each 
independent direction of the phase space, so the chaotic 
properties of many-body systems can be characterized by 
an ordered set of Lyapunov exponents, the so-called Lya- 
punov spectrum {A' 1 ', A^ 2 -*, • • • } where A'- 1 ) > A'- 2 -' > • •• . 
The Lyapunov spectrum is connected to the contraction 
rate of the phase space volume (roughly speaking, the 
dissipation rate) through the sum of all the Lyapunov ex- 
ponents, allowing the calculation of transport coefficients 
from the Lyapunov spectra (HQ- El The conjugate pair- 
ing rule for Lyapunov spectra reduces the calculation of 
the sum of all the Lyapunov exponents to the sum of just 
one pair of Lyapunov exponents 0, IE IS ■ The set of 
all positive Lyapunov exponents specifies the natural in- 
variant measure |j| @ > which is used to calculate various 
quantities using periodic orbit theory [lCfl and led to the 
first form of the fluctuation theorem |lllll2l fl3| . On the 
other hand, recently, much attention has been paid to 
individual Lyapunov exponents and their Lyapunov vec- 
tors for many-body systems. As each Lyapunov exponent 
has the dimensions of inverse time, the Lyapunov spec- 
trum can be regarded as a time-scale spectrum. From 
this point of view, Lyapunov exponents with small abso- 
lute values are connected to large and macroscopic time- 
scale behavior of many-body systems, and in this region 
the wave-like structure of Lyapunov vectors, known as 
the Lyapunov modes, is observed 0, 0, 0, 0, Il8| . 
The Lyapunov mode is a reflection of a collective move- 
ment (phonon mode) of many-body systems, and comes 
from dynamical conservation laws and translational in- 
variancc 

EES HI El On the other hand, large Lya- 
punov exponents are dominated by small and microscopic 
time-scale movement, and in this region the spatially lo- 



calized behavior of Lyapunov vector, the so-called Lya - 
punov localization, appears pi El iM El El El |2S| . 
For the largest Lyapunov exponent of many-body chaotic 
systems , a nal y tical calculations have been attempted 
[2I I3I I3H El l33| . The time-scale separation in many- 
body systems is crucial to extract a macroscopic dynam- 
ics from microscopic many-body dynamics, and the Lya- 
punov spectrum allows us to discuss it dynamically. 

As one of the features of Lyapunov vectors for many- 
body systems, the Lyapunov localization appears as a 
behavior in which the Lyapunov vector components for a 
few particles are significantly larger than the other com- 
ponents. Moreover, the localized region of a Lyapunov 
vector moves as a function of time. The magnitude of 
the localization of each Lyapunov vector can be mea- 
sured quantitatively by the Lyapunov localization spec- 
trum, which is defined as a set of exponential functions 
of entropy-like quantities for the normalized amplitudes 
of the Lyapunov vector components [23, OS |3a] . We have 
previously reported that the Lyapunov localization spec- 
tra show a bending behavior at low density, and its con- 
nection with kinetic theory properties (e.g. the Krylov 
relation for the largest Lyapunov exponent, and the mean 
free time being inversely- proportional to density, etc.) 
is discussed for many-hard-disk systems [28|, |34| . How- 
ever, the Lyapunov localization spectrum requires taking 
a time-average and characterizes only the static localiza- 
tion of the Lyapunov vectors. The physical meaning of 
the movement of the localized regions of Lyapunov vec- 
tors is not clearly understood. 

The principal aim of this paper is to discuss dynam- 
ically the movement of the localized region of the Lya- 
punov vectors in many-hard-disk systems at low density. 
To discuss this problem we use the fact shown in Ref. 
|28| that at low density only two particle components of 
the normalized Lyapunov vector for large Lyapunov ex- 
ponents have a non-zero value. This is due to the short- 
range of particle interactions in a many-hard-disk system. 
The movement of the localized region of the Lyapunov 
vector appears to be a series of jumps or hops, so we can 
introduce a hopping rate to describe the dynamics. To 
simplify the problem, in this paper we consider quasi- 
one-dimensional systems, in which the system width is 
so narrow that disks always remain in the same order 
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EEl IH HI We show that this hopping rate depends 
on the hopping distance, and is a decreasing function of 
the hopping distance. This implies that there is a spatial 
correlation among localized regions of Lyapunov vectors. 

We explain the hopping-distance dependence of the 
hopping rate in many-hard-disk systems in two ways. In 
the first approach we use a simple model expressed as an 
accumulation of bricks. Here, the hopping of the localized 
region of the Lyapunov vectors is expressed as a change 
in the position of the highest brick site. This model is 
a one-dimensional version of the so-called clock model, 
which has been used to calculate Lyapunov exponents 
for many-hard-disk systems [3lL 13^ . 133 ] . In this paper we 
demonstrate that this model can reproduce the largest 
Lyapunov exponent for quasi-one-dimensional hard-disk 
systems. As the second approach to the hopping behavior 
of the localized Lyapunov vectors, we propose an analyt- 
ical method to calculate the hopping rate from the sum 
of probability distributions for the brick height configu- 
rations between two separated highest brick sites. Using 
this analytical approach we can also discuss the relation 
between the hopping rate and the Lyapunov exponent. 
Hopping rates calculated by these two approaches are in 
good agreement with the ones for quasi-one-dimensional 
hard-disk systems. 

The outline of this paper is as follows. In Sec. [H] 
we introduce the quasi-one-dimensional hard-disk sys- 
tem, and show the localized behavior of the Lyapunov 
vector corresponding to the largest Lyapunov exponent 
at low density. In Sec. IIHI we introduce the hopping 
rate of the localized region of the Lyapunov vectors, and 
show the hopping-distance dependence for many-hard- 
disk systems. In Sec. IIVI we discuss the hopping rate 
using a brick accumulation model. In Sec. |V|we propose 
an analytical expression for the hopping rate. Section 
IVII is our conclusion and some remarks. In Appendix 
1X1 we discuss the technical details of the calculation of 
the hopping rate. In Appendix 151 we give a microscopic 
derivation of the brick accumulation model. 



II. QUASI-ONE-DIMENSIONAL SYSTEM AND 
LOCALIZATION OF LYAPUNOV VECTORS 

The system which we consider in this paper is a quasi- 
one-dimensional system consisting of N hard-disks in pe- 
riodic boundary conditions. All of the particles are iden- 
tical with radius R and mass M, and the shape of the sys- 
tem is rectangular with the length L x and the width L y 
satisfying the inequality 2R < L y < AR. The schematic 
illustration of the system is given in Fig. ^ in which 
we number particles 1, 2, • • • , N from the left to right in 
this system. For the actual numerical results shown in 
this paper, we used: the radius of a particle R — 1, the 
mass of a particle M — 1, the total energy of the sys- 
tem E = N, the system width L y = 2R(1 + 1CT 6 ), and 
the system length L x — NL y (l + d) with the constant d 
controlling the density p = NirR 2 /(L x L y ). In the quasi- 
one-dimensional system, the particle interactions are re- 
stricted to nearest-neighbor particles only, so particles 
remain in the same order. These features require less cal- 
culation effort and a simpler representation of results for 



FIG. 1: A quasi-one-dimensional system consisting of hard- 
disks with radius R. The length L x and the width L y of 
the system control the density. The width L y satisfies the 
inequality 2R < L y < 47?, so that the disks remain in the 
same order and particles can be numbered 1, 2, ■ • ■ , N from 
the left to right. 



quasi-one-dimensional system compared with fully two- 
dimensional systems. The quasi-one-dimensional system 
has already been used to investigate the localized behav- 
ior of Lyapunov vectors I2SL 1.341 . the wave-like structure 
of Lyapunov vectors Hil |21| and the transition be- 
tween quasi-one-dimensional and fully two-dimensional 
systems [3fij . 

The dynamics of Lyapunov vectors in many-hard-disk 
systems is well established, and readers should refer to 
the references, for example, Ref. |4ll |. for more detailed 
discussions. In many-hard-disk systems the dynamics is 
separated into a free-flight part and a collision part, and 
the free-flight part of the dynamics is integrable. This 
property allow us to express the dynamical evolution as 
a simple multiplication of time-evolutional matrices for 
the free-flight dynamics and the collision dynamics, lead- 
ing to a fast and more accurate numerical simulation than 
for soft-core interaction models. For numerical calcula- 
tions of the Lyapunov vectors shown in this paper, we 
use the algorithm developed by G. Benettin, et al [37} 
and I. Shimada and T. Nagashima (also see Refs. 
[3^. lifTj). This algorithm is characterized by intermit- 
tent (e.g. after every collision) re-orthogonalization and 
renormalization of Lyapunov vectors, preventing a diver- 
gence of the amplitude of the Lyapunov vectors. 

In this paper we use the notation <5T ,( "' ) (t) = 
OSrjt^CO.&r^i),-. • ,5r<g\t)) for the Lyapunov vec- 
tor corresponding to the n-th Lyapunov exponent \( n ' at 
time t. Here, ST^ (t) is the Lyapunov vector component 
contributed by the j-th particle in the n-th Lyapunov ex- 
ponent at time t. To express the localized behavior of the 
Lyapunov vectors we introduce the quantity 7^™ (t) as 

Mf>(t) 2 

E «f >(*) 
fc=i 

which is the normalized amplitude of the Lyapunov vec- 
tor component for the j-th particle for the n-th Lya- 
punov exponent at time t. The localized behav- 
ior of Lyapunov vector, namely the Lyapunov localiza- 
tion, is the phenomenon where only a few of the 7^™' (£) , 
j = 1, 2, • • ■ , N have a non-zero value at any time t. 

Fig. Elis an example of the Lyapunov localization. It is 
a graph of the normalized amplitude 7^ of the Lyapunov 
vector components for the j-th particle corresponding 
to the largest Lyapunov exponent AW as a function of 
the particle index j and the collision number n t (~ t/r 
with time t and the mean free time r). The system is 
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FIG. 2: The normalized amplitude 7^ of the Lyapunov vec- 
tor particle component corresponding to the largest Lyapunov 
exponent as a function of the collision number n t and the par- 
ticle index j in a quasi-one-dimensional hard-disk system with 
d — 10° and TV = 50. On the base of this graph is a contour 
plot of 7^' at the level 0.2. 



FIG. 3: The normalized hopping rates Pjv(/i)/Pjv(l) for 
quasi-one-dimensional systems of different sizes at the same 
density (d = 10 3 ) as a function of the hopping distance \h\. 
The number of particles is TV = 25 (circles), 50 (triangles), 
and 75 (squares). The graphs are log-log plots with error 
bars are given by |Pjv(/i) — Pjv(— ft)|/Pjv(l). The lines are 
fits of the numerical data to functions based on Eq. (|5J with 
power functions as the asymptotic form of the hopping rate. 



a quasi-one-dimensional system consisting of 50 hard- 
disks and with d = 10 5 (the density p » 7.85 x lfT 6 ). 
From this figure we recognize that the non-zero compo- 
nents of the Lyapunov vector are concentrated at two 
nearest-neighbor particles. This characteristic is shown 
for hard-disk systems at low density and for large Lya- 
punov exponents in general, using the Lyapunov localiza- 
tion spectrum [2£||34|]- Moreover, we observe that such 
localized regions of Lyapunov vector components move 
with time in discrete jumps or hops, in Fig. [21 This 
characteristic hopping has already been shown in Ref. 
[28| and that such spatial hopping of the non-zero com- 
ponents of 7; (t) are caused by some particle-particle 
collisions. However, not every collision causes a hopping 
of the localized Lyapunov vector. In this sense, parti- 
cle collisions themselves are not sufficient to explain the 
hopping movement of the Lyapunov localization. 

The Lyapunov localization appears in the Lyapunov 
vectors corresponding to large Lyapunov exponents in 
general, but in this paper for simplicity we consider only 
the Lyapunov vector ST^ corresponding to the largest 
Lyapunov exponent . 



III. HOPPING RATE OF LOCALIZED 
LYAPUNOV VECTORS 

An advantage of the quasi-one-dimensional system in 
investigations of Lyapunov localization is that the move- 
ment of particles in the transverse direction are sup- 
pressed, and roughly speaking, the particle sequence cor- 
responds to the particle's position. Noting this feature, 
in this section we describe the hopping behavior of spa- 
tially localized Lyapunov vectors as the hopping of par- 
ticle indices whose defined by Eq. (JTJ have non-zero 
values. 

As shown in Refs. [28l l3~3 | the particle indices with 
non-zero 7^ are a pair of nearest-neighbor particles in 
the low density limit, and we can introduce the hopping 
distance h at each collision as the change of particle in- 
dices. In this definition the hopping distance h is an 
integer number satisfying the inequality —[N/2] < h < 
[N/2], where [x] is the integer part of the real number 
x. It should also be noted that we take the particle 
index j as the one equivalent to the index j ± N , be- 
cause of periodic boundary conditions, so the hopping 
distance h± N is equivalent to h. More technical details 
of the calculation of the hopping distance h are given 
in Appendix [S] Using this hopping distance h, in nu- 
merical simulations, we count the number TVr(/i) of hops 



with hopping distance h in a time-interval T, and we de- 
fine the hopping rate P^{h) a $ a function proportional 
to Nx(h) as T — > 00. In this paper we use the hop- 
ping rate P N {h)/P N (l) = limr^oo N T {h)/N T {l) nor- 
malized by Pjv(l). From the reflection symmetry of the 
quasi-one-dimensional system in the longitudinal direc- 
tion, the hopping rate Pjv(/i) must be symmetric, namely 
Pn(It-) = Pff(—h). We use this hopping rate to quanti- 
tatively discuss the hopping dynamics for localized Lya- 
punov vectors. 

Figure |3] shows log-log plots of the hopping rates 
P/v(/i)/P/v(l) normalized by P/v(l) as a function of \h\ in 
a quasi-one-dimensional hard-disk system with d = 10 3 
at different numbers of particles; N = 25 (circles), 50 (tri- 
angles), and 75 (squares). Noting the symmetric prop- 
erty P N (-h) = P N {h), we use \P N (h) - P N (-h)\/P N (l) 
as error bars in Fig. |21 It is clear from Fig. |21 that 
the hopping rate Pjv(|/i|) decreases as \h\ increases. This 
implies that there is a spatial correlation among the local- 
ized regions of a Lyapunov vector, rather than a random 
hopping. 

The turn-up in the tail of the normalized hopping 
rate Pjv(/i)/Pjv(1) m Fig- 01 can be explained as an 
effect of periodic boundary conditions. Under periodic 
boundary conditions the hopping distance h + jN for 
any integer j is observed as the hopping distance h in 
the quasi-one-dimensional system consisting of TV hard 
disks. Therefore, using the hopping rate Poo(h) for the 
thermodynamics limit (N — > 00 at a fixed density) for 
h = 0, ±1, ±2, ••• ± 00, the hopping rate P/v(/i) for a 
finite system should be represented as 

00 

Pn(H)= P oo{h + jN) (2) 

j=— 00 

for/i= -[iV/2],-[TV/2] + l,--- ,[N/2]. The terms on the 
right-hand side of Eq. (J2J with j ^= cause the turn-up 
in the tail of the hopping rate P/v(/i) for a finite system. 
In this sense we can explain the TV-dependence of the 
hopping rate Pj^(h) using an TV-independent asymptotic 
distribution Poo(h). To make this explanation convincing 
we fitted the numerical data for Pjv(/i)/P/v(1) shown in 
Fig. CDto the function f(h) = a[h +J^=i (\h+jNf+\h- 
jN\P)] assuming that the decay of Poo(h) = a\h\P with 
fitting parameters a and /?, and neglecting the higher 
order small terms for \j\ > 3. Here, the fitting lines are 
dotted for TV = 25, the broken for TV = 50 and solid for 
N = 75. Note that these sets (a,/3) of fitting parameter 
values, summarized in Table [I] are almost independent 
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TABLE I: Fitting parameters for the power function (h) = 
a\Kf for quasi-one dimensional systems of different numbers 
of particles N with d — 10 3 . Notice that the coefficients a 
and are essentially independent of N. 



N 


a 


P 


25 


0.718 


1.67 


50 


0.788 


1.72 


75 


0.753 


1.71 



FIG. 4: The brick accumulation model expressing the dy- 
namics of the amplitude of the Lyapunov vector component 
for each particle in a quasi-one-dimensional system. 



of the number of particles N. As shown in Fig. these 
fitting lines nicely reproduce the values of the numerical 
data, especially the turn-ups. 

We can also calculate the hopping rate Pjsr(h) at dif- 
ferent densities, as long as the density is low enough 
so that only a few of the normalized Lyapunov vec- 
tor component amplitudes {7j- }j have non-zero val- 
ues. In Appendix ^ we show that the normalized hop- 
ping rate P/v(/i)/P/v(l) is almost density-independent at 
least for 10 3 < d < 10 5 , (namely in the density region 
7.85 x 1CT 6 < p < 7.85 x 1CT 4 ). However, a subtle 
increase of the normalized hopping rate as the density 
decreases is recognizable in this low density region. 



IV. BRICK ACCUMULATION MODEL 



As shown in Sec. IHII the hopping rate of the localized 
region of the Lyapunov vector is correlated spatially. In 
this section we explain this characteristic using a simple 
one-dimensional model, which we call the brick accumu- 
lation model. 

A schematic illustration of the brick accumulation 
model, or simply the brick model, is given in Fig. It is 
a one-dimensional model with N sites in the horizontal 
direction. The bricks are dropped at random and occupy 
a pair of neighboring sites. Each brick has width two and 
the height 1 and they accumulate at sites without over- 
lap. The dynamics of the brick model is described using 
the brick height JCj (n) at the j-th site after n bricks have 
been dropped. The total brick height JCj (n) takes on in- 
teger values and its dynamics is expressed as follows. If 
the n-th brick is dropped on sites j n and j n + 1 then 



K.i{n) 



max{/Cj 7l (n — l),/Cj„+i(n — 1)} + 1 
if I € {j„, j„ + 1} 

Ki{n-1) if l<£{j n ,j n + l} 



(3) 



(noting that the particle index N + 1 is equivalent to 1). 
Here the site number j n is chosen randomly on [1, N] for 
each n. 



One may notice that the brick accumulation model is 
a one-dimensional version of the so-called clock model 
[3TI 135. l33| , which has been used to calculate the Lya- 
punov exponents for many- hard-disk systems. In the 
clock model, the brick height JCj(n) is usually called the 
clock value of the j-th particle after the n-th particle- 
particle collision, and the dynamics Q is expressed as 
an adjustment of the clock values of colliding particles. 
However we use the name brick accumulation model for 
the one-dimensional version of the clock model, because 
in a simple image of brick accumulations, it is easily to 
visualize the configuration of brick heights which is es- 
sential for the analytical approach to the hopping rate 
discussed in the next section^ although this image is ap- 
plicable only for the one-dimensional case. The image of 
brick accumulations also helps us to easily recognize the 
similarity of this model to ballistic aggregation models, 
whose scaling properties have been studied analytically 
and numerically 0> 03 ■ 

In the brick model described by the brick height JCj (n) , 
the site number j corresponds to the particle index, and 
the number n of dropped bricks corresponds to the col- 
lision number n t for the quasi-one-dimensional hard-disk 
system. Moreover, the brick height ICj(n) itself is con- 
nected to the Lyapunov vector component amplitude 



\srf\t)\ by 



5rf\0) cxp{-/C 3 Mlnp} (4) 



for the particle index j and the collision number nt with 
the mean free time r in the asymptotic limit of low den- 
sity p — > 0. The derivation of the relation (0} from mi- 
croscopic dynamics for hard-particle systems is given in 
Appendix [5] 

From the relation J2J| between the brick height and 
the Lyapunov vector, the amplitude of the localized Lya- 
punov vector components with the highest brick height 
have much larger values than the other components be- 
cause of the huge factor — In p appearing in the exponent 
of Eq. 0} at low density p « 1. Therefore, to a good 

in) 

approximation, the normalized amplitude given by 
Eq. has non-zero components only for particles cor- 
responding to these largest amplitudes, as shown in Fig. 
12 The combination of Eqs. © and (@J also explains 
why the amplitudes of the localized components of the 
Lyapunov vectors for nearest-neighbor particles take al- 
most the same value, thus appearing as a flat top in the 
localized regions in Fig. [2J After all, the localized re- 
gion of the Lyapunov vector represented in Eq. 0} is 
given as the site indices whose brick height JCjijit) take 
the maximum value after the n-th brick is dropped, and 
such sites with the highest brick height can be calculated 
using the brick model dynamics © only, without refer- 
ring to the hard-disk dynamics. Based on these features 
we can calculate the hopping distance h for the brick ac- 
cumulation model (see the end of Appendix lAl for more 
detail of the hopping distance for the brick model), and 
therefore the hopping rate P/v(/i), similarly to the quasi- 
one-dimensional hard-disk system. 

Figure [S] shows the normalized hopping rates 
P/v(/i)/Pat(1) for the brick model with 50 sites (cir- 
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FIG. 5: Normalized hopping rates Pjv(ft.)/Pzv(l) for the brick 
model (circles) and the quasi-one-dimensional system with 
d — 10 5 (squares) as a function of the absolute value of hop- 
ping distance \h\ for TV = 50 as log-log plots. The error bars 
are given by |Pjv(/i) — Pn(— ft)|/Piv(l). The broken line is the 
analytical expression for the hopping rate discussed in Sec. 

El 



FIG. 6: The sum ^2j =1 ICj(n) of brick heights /C,(n) as a 
function of the number n of dropped bricks in the brick accu- 
mulation model with N = 50 in the n interval [0,20979999]. 
Inset: The same graph but enlarged in the much smaller n 
interval [1000, 1050]. The broken line (which is almost on the 
data points of the main graph) is a fit of numerical data to a 
linear function. 



TABLE II: Particle number (N) dependences of the ratio 
Pjv(0)/Pjv(l) of the hopping rates at the hopping distances 
h = and 1 in the quasi-one-dimensional hard-disk systems 
with d — 10 3 and in the brick accumulation model. 





Hard-disk model 


Brick model 


N 


Pv(0)/Pv(l) 


Pv(0)/Pv(l) 


25 


13.5 


17.1 


50 


26.8 


37.9 


75 


39.6 


58.8 


100 


52.1 


79.7 



cles) and for a quasi-one-dimensional hard-disk system 
consisting of 50 hard-disks and d = 10 5 (squares). 
Agreement between the brick model and the quasi-one- 
dimensional hard-disk system for the normalized hopping 
rate Pjv(/i)/-Pjv(l) f° r \h\ < 1 is satisfactory. Numerical 
simulations of quasi-one-dimensional hard-disk systems 
show that the normalized hopping rate P/v(ft)/P/v(l) in- 
creases very slightly as the density decreases as shown in 
Appendix lAl so the deviation of P/v(ft)/P/v(l) in tails in 
the brick model and the hard-disk system in Fig. [S]may 
be a finite density effect. 

So far we have discussed the hopping rate Pat (ft.) for 
non-zero hopping distances ft ^ 0. Now we discuss the 
hopping rate Pjv(0) for zero hopping distance ft — 0. Ta- 
blelTTI shows the ratio P/v(0)/P/v(l) for the hopping rates 
at ft = and 1 in the quasi-one-dimensional hard-disk 
system with d = 10 3 and the brick accumulation model 
for N = 25, 50, 75 and 100. As shown in this table, the ra- 
tio Pjv(0)/Pjv(1) depends on the number of particles N, 
and roughly speaking it is proportional to N. To explain 
this linear dependence of Pjv(0)/P/v(1) with respect to 
N, we note that in the brick model the hopping distance 
1 occurs only when a new brick is dropped at one (and not 
both) of the sites with the highest brick height, meaning 
that the hopping rate Pjv(1) should be approximately 
inversely proportional to N . On the other hand, the 
hopping distance occurs when a new brick is dropped 
at a site which does not have the highest nor the sec- 
ond highest brick height, or is dropped at both the sites 
which are currently the site of the highest brick height, 
meaning that the hopping rate Pjv(0) should be, roughly 
speaking, independent of N. These considerations for 
Pjv(0) and P/v(l) explain why the ratio P/v(0)/P/v(l) is 
proportional to N. However we need a more detailed 
consideration to explain the difference between the coef- 
ficient k in the relation Pjv(0)/Pjv(1) ~ kN between the 



quasi-one-dimensional hard-disk system and the brick ac- 
cumulation model. 

Before finishing this section, we discuss one more prop- 
erty of the brick accumulation model, which we will use 
in the next section. Fig. HJ1 is the sum Y^j=i^-j{ n ) °f 
brick heights as a function of n. Note that this sum must 
increase at a speed of more than 2 per dropped brick, 
because accumulated bricks capture spaces below them 
which cannot be occupied by further dropped bricks. It is 
clear from this figure that this sum increases linearly, and 
a fit of the data to the function Y2j=i ( n ) ~ ot+fin with 
fitting parameters a and (3 leads to the values a ~ 0.00 
and (3 « 3.99 ^?|- (Note that data points in the main 
graph of Fig. look exactly like this fit because of the 
large scale. In the inset to Fig. H^we show the graph of 
£j=i fcj( n ) as a function of n on a much smaller scale, 
to show its fluctuating behavior.) Therefore, on average, 
each dropped brick adds a contribution of 4 to the sum 
53j=i meaning that each brick occupies not only 

its own 2 spaces but also captures 2 empty spaces below 
it. 

It is important to note that from Eq. @ the brick 
height !Cj(n) dominates the exponential growth rate of 
the Lyapunov vector component amplitude. Using this 
feature of the brick heights, their sum YljLi &j ( n ) i s con ~ 
nected to the largest Lyapunov exponent A' 1 ' as 



A«~- lim — 



1 N 



In p 



(5) 



with the mean free time r and the density p in the asymp- 
totic limit of low density. [More detailed discussion for 
the formula JSJ is given in Appendix [5]] As a numerical 
check of the formula JSJ) we show in Fig. the largest 
Lyapunov exponent A*- 1 -* as a function of the density p in a 
quasi-one-dimensional hard-disk systems and in the brick 
accumulation model using the formula (JSJ with N = 50. 
Here, to calculate the largest Lyapunov exponent from 
Eq. (J5J we used the relation J2j=i ^-j( n ) ~ 4n as shown 
in Fig. and the values of the mean free time r and 
the density p of the quasi-one-dimensional system whose 
Lyapunov exponents are plotted in Fig. [71 and the data 
points are connected by a dashed line for ease of visibility. 
Figure [7J shows that Eq. (J5J reproduces successfully the 
values of the largest Lyapunov exponent for the quasi- 
one-dimensional hard-disk systems, not only in the limit 
of low density but also at relatively high density such as 
p < 0.3. It may be noted that a linear dependence of 
the sum J2 7=1 ICj (n) of brick heights with respect to n 



6 



FIG. 7: The largest Lyapunov exponent A' 1 ' as a function 
of density p in a quasi-one-dimensional hard-disk system with 
N = 50 as a log-log plot. The error bars are given by |A^ — 
\( 4N ) | which must be zero by the conjugate pairing rule for 
Hamiltonian systems. The broken line is the largest Lyapunov 
exponent given by the brick accumulation model with N = 50. 



is necessary to get a finite value of the largest Lyapunov 
exponent A' 1 ) by Eq. (J5J. 



FIG. 8: Schematic illustration for a brick-height config- 
uration (similar to Fig. [I] for the brick model) connect- 
ing two highest sites (particle indices) [i and p: + h with 
hopping distance h in the case of h > 2. Here, —ki 
is the brick-height difference of the sites fj, + I — 1 and 
(i + I (I = 1, 2, • • • ,h — 1). The probability distribution for 
the configuration with ki, I — 1,2, ••• ,h — 1 is given by 
A(-fci)A(-fc 2 )A(-fc 3 ) ■ ■ ■ A(-fc h _ 1 )A(J]^r i 1 kj) with the dis- 
tribution A(k) of brick-height differences k between nearest- 
neighbor sites. The hopping rate Pao(h) is calculated as the 
summation of this probability distribution over possible val- 
ues of ki, I = 1, 2, • • ■ , h — 1. 



V. ANALYTICAL EXPRESSION FOR THE 
HOPPING RATE 

In this section we discuss another approach to the hop- 
ping rate of the localized region of the Lyapunov vectors. 
This approach is inspired by some of the characteristics 
of the brick accumulation model discussed in Sec. IIVI al- 
though we greatly simplified the brick model by omitting 
some other aspects of the model. The advantage of this 
approach is that we can get an analytical expression for 
the hopping rate, while to a good approximation it still 
reproduces the hopping rate for hard-disk systems. An- 
other important point in this approach is that it connects 
the dynamics of the localized region of the Lyapunov vec- 
tors with a static property, the probability distribution of 
brick-height differences between nearest-neighbor sites. 

In the brick accumulation model, a hop of the highest 
brick site occurs when two (non-nearest-neighbor) sites 
have the same brick height. Noting this characteristic 
we consider the probability distribution Poo(h) under the 
constraint that the two sites /j, and n + ft have the highest 
brick height. We require that there is no other highest 
site between the two highest sites /i and fi+h for \h\ > 2; 



Ki{n) < KJn) 



[= /C M+7l (n)], 

+ 1, (i + 2, • • • , n + ft — 1 
for ft > 2 
/i + ft + 1 , /i + ft + 3, • • • ,fJ,— l 
for ft < -2, 



(6) 



Using this probability distribution Poo (ft), we can es- 
timate the hopping rate P>o(ft) in the thermodynamic 
limit as the one proportional to Poo (ft): Poo (ft) oc Poo(h). 
Then we can calculate the hopping rate P/v (ft) for a finite 
size system using Eq. P|l. apart from a constant factor. 

To calculate the probability distribution Poo (ft) we in- 
troduce the distribution A(fc) of brick-height differences 
k between nearest-neighbor sites. These two distribution 
functions are connected by 

£»w = E E ••• E 

k% ko, fc|h|-i 

(fei>l) (fei+fe 2 >l) fej>l) 



'\h\-\ 

:A(-fc 1 )A(-fc 2 )-.-A(-fc| fthl )A ( k i 

3=1 



In Eq. (JJJ, — ki is a brick- height difference of near- 
est neighbor sites (/ = 1, 2, • • • , \h\ — 1). The prob- 
ability distribution for the specific brick-height con- 
figuration with ki, I = 1, 2, • • • , \h\ — 1 is given by 
A(-h)A(-k 2 )A(-k 3 ) ■ ■ ■ A(~fc|^|_ 1 )A(E;"Jr 1 and 
the hopping rate Poo(h) is calculated as the summation 
of this probability distribution over possible values of ki , 
/ = 1,2, ••• ,h — 1, leading to Eq. Q. Figure El is a 
schematic illustration of a brick- height configuration con- 
necting two highest sites and the brick-height differences 
—kj of nearest neighbor sites in the case of h > 2. In Eq. 
(0) , the lower bounds of the sums on the right-hand side 
of Eq. J7J come from the condition |J?jJ. Note that from 
Eq. JJJ the distribution Poo(h) is an even function of h, 
namely Poo(-/i) = AoW- 

Now, for simplicity, we assume that the distribution 
function A(k) can be expressed as an exponential func- 
tion 



A(fc) = Wcxp{-? ? |fc|} 



(8) 



with constants W and rj (> 0). (We will discuss the 
validity of this assumption later.) Inserting Eq. |JSJ into 
Eq. 0, and replacing the sums over ki in Eq. Q with 
the ones over k[ = ki + 1, we obtain 



k=0 



for | ft | 



Poo (ft) = < 



-f-oo -t-oo -t-oo 

ki—O k°2 = — k\ k^ — — k\ — k2 

-j-oo 

|M-2 

fc |/i|-l = _ E fc 3 
3=1 

for |ft| > 3 



where T is defined by 

T = cxp{— 2rj} 



(9) 



(10) 



(7) 



and 0(x) is the Heaviside function taking the value 1 
for x > and the value for x < 0. The summations 
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appearing in Eq. JHJ can be carried out successively, and 
we obtain 



yv\ h \T 

f °°W= (1 _ T) | h |-i "(N) 

where the function 0(/c) of k is can be written as 
0(2) = 1, 

n(3) = i + T, 

0(4) = 1 + 3T + T 2 , 

0(5) = (1 + T)(1 + 5T + T 2 ), 



0(6) 



1 + 10T + 20T 2 + 10T J + T 



(11) 



(12) 
(13) 
(14) 
(15) 
(16) 



0(7) = (1 + T)(1 + 14T + 36T 2 + 14T 3 + T 4 ), (17) 

and so on. Eq. withEqs. Q], flU, 0, fTB) . 

(|17|l . etc., gives an analytical expression for the hopping 
rate, for example, using the relation Poo(M/P=o(l) = 
(1). 

The coefficients W and appearing in Eq. (JSJ) can be 
determined from the two sum rules: 



E A ( fc ) = 

A:— — og 

E i fc i A ( fc ) = A/c > 



(18) 
(19) 



k=— oo 



where AAC is the mean value of the absolute value of the 
brick height difference between nearest-neighbor sites. 
The first condition H18|) is the normalization of the prob- 
ability distribution A(fc) which leads to 



W 



exp{r?} - 1 



exp{^} + 1 

Using Eq. (|2U[1 the second condition (|19|l gives 



r) = In 



1 + V1 + A/C 2 
A/C 



(20) 



(21) 



satisfying the inequality r\ > 0. Inserting Eq. (|21|l into 
Eqs. l(Tfj|) and (|2"0j) we obtain 



T 



A/C 



1 + VI + A/C 2 



1 + VI + A/C 2 - A/C 
1 + VI + A/C 2 + A/C ' 



(22) 
(23) 



From Eqs. IL'l't and l|23|l . there is only one parameter 
A/C remaining to determine the hopping rate using Eq. 

To estimate the value of A/C, we use a property of 
the brick accumulation model discussed at the end of 
Sec. IIVI Previously, we showed an approximate re- 
lation ~Y^=ifcj{n) ~ 4n for the sum of brick heights, 
which means that in the brick accumulation model each 
dropped brick gives a mean contribution of 4 to this sum. 
This contribution consists of a contribution of 2 as the 
space occupied by a brick itself and another 2 as empty 



space below the brick which can now not be occupied by 
other bricks. This implies that the averaged brick-height 
difference between nearest-neighbor sites is about 2, so 
that 



A/C 



(24) 



We use this value to calculate the hopping rate based 
on Eq. (|llfl . One may notice that from the relation 
J2jLi ^3 ( n ) ~ (A/C + 2)n and the formula we obtain 



ak ~ *r a« 

lnp 



(25) 



in the limit of low density. From the relation l|25|l between 
the parameter A/C specifying the hopping rate and the 
largest Lyapunov exponent A^ , the hopping rate of the 
localized region of the Lyapunov vectors is connected to 
the largest Lyapunov exponent. 

Before comparing the hopping rate based on the an- 
alytical expression l|ll|) with the ones for hard-disk sys- 
tems and the brick accumulation model, we note that 
this analytical expression for the hopping rate may not 
be appropriate for small hopping distances \h\. In partic- 
ular, it does not give the correct value of P/v(l), because 
in the brick model the hopping distance h = ±1 does not 
appear from separated highest sites with the same brick 
height, the assumption used to derive Eq. (|TT)l . There- 
fore it is not appropriate to calculate the hopping rate 
normalized by Pat(1) from this analytical expression and 
to compare it with the numerical results. The hopping 
rate Pat (±2) from Eq. i|ll|) may also be problematic, be- 
cause in the brick model the hopping distance h = ±2 
occurs when 4 consecutive sites have the highest brick 
height, while to derive Eq. (|llf> we assumed that h = ±2 
occurs when non-consecutive separate sites /i and fi + h 
have the same highest brick height. Based on these con- 
siderations we do not calculate the value P/v(±l) from 
the analytical approach in this section, and plot the hop- 
ping rate so that Pat (5) by Eq. (fTT|l gives the same value 
as that from the brick model. 

In Fig. [3] for the normalized hopping rate 
P N (h)/P N (l), we plotted the function = 
Pn (h)P N {5)/[P N (5)P N (l)} using the value P N {5)/P N (1) 
of the brick model with P N (h) = P^/i) + Poa(N - h) 
using the analytical expression Poo(h) given by Eq. l(TT|) 
for \h\ > 2. Note *(5) = Pv(5)/Pv(l) so that co- 
incides with Pat(/i) / Pat(1) of the brick model at \h\ = 5. 
Here, the value of the hopping rate values are given at 
integer values of h, but we connect them with a broken 
line for ease of visibility. Figure shows that the analyt- 
ical expression (lllfl for the hopping rate reproduces the 
hopping rate for the brick accumulation model as well 
as the quasi-one-dimensional hard-disk system to a good 
approximation. It may be noted that for this plot we 
used the first two dominant terms on the right-hand side 
of Eq. (J2J), and part of the small deviation of the hopping 
rate in the tail between the brick model and the analyt- 
ical expression should come from the omission of higher 
order terms in Eq. J2J. 

We notice that the approach in this section is sim- 
ple enough to get an analytical expression for the hop- 
ping rate but it is not completely consistent with the 
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FIG. 9: The distribution A(k) of the brick-height differences 
k between nearest-neighbor sites as a function of \k\ in the 
brick accumulation model with iV = 50 as a linear- log plot. 
Here, A(fc) is normalized by X]fc=-2iv A-(fc) = The distribu- 
tion A(fe) is an even function of k, and the error bars in this 
graph are given by \A(k) — A(— k)\. The line is an exponen- 
tial function used to obtain the analytical expression for the 
hopping rate given by Eqs. JHJ, 1221 . and t3H , 



FIG. 10: The distribution D(a ■ Ap) of the quantity a ■ Ap 
with the normalized collision vector cr and the momentum 
difference Ap of colliding particles just before the collision; 
for the general case (dashed line) and for the case in which 
non-zero hopping of a localized Lyapunov vector occurs. The 
system is a quasi-one-dimensional system with d — 10 3 and 
N = 25. 



brick accumulation model discussed in the previous sec- 
tion IIVI Previously, we have mentioned the irrelevance 
of Eq. as a description of a small hopping rate in 

the brick model. Actually, in the approach of this sec- 
tion we omitted the essential characteristic of the bricks 
as components of the brick accumulation model, except 
for the property l|24|) , and treated the model components 
as blocks (or half bricks). As another example, we show 
in Fig. EJthe numerical result for the distribution A(£;) of 
brick-height differences k between nearest-neighbor sites 
in the brick accumulation model with N = 50. Here, 
A(fc) is normalized by Y^t=-2N = 1> instead of Eq. 
(|18fl . because we cannot calculate A(fc) in |fc| — > +oo 
numerically. In this figure we added the exponential dis- 
tribution © using Eqs. (HJ}, H3 and Figure H 
shows that the probability distribution A(fc) does not co- 
incide with an exponential distribution JHJ, although it 
may be justified as a first approximation. (On the other 
hand, the numerical evaluation of AK, from the distribu- 
tion A(fc) as a numerical result in Fig. gives the value 
1.99, then Eq. I|24|l is still justified.) On another point, 
in the brick model the highest brick sites appear as a 
pair of nearest-neighbor sites, but we did not take into 
account this characteristic in the analytical approach in 
this section. Despite the omission of characteristics of the 
brick accumulation model, the analytical approach in this 
section reproduces the hopping rate for many-hard-disk 
systems reasonably well, and it suggests that this ap- 
proach still keeps enough of the essential characteristics 
that describe the dynamics of the Lyapunov localization. 



VI. CONCLUSION AND REMARKS 

In this paper we discussed the dynamics of the spa- 
tially localized region of the Lyapunov vector correspond- 
ing to the largest Lyapunov exponent in a quasi-one- 
dimensional hard-disk system. To discuss the dynamics 
of the localized region of the Lyapunov vector we intro- 
duced a hopping rate for the localized region, and showed 
that the hopping rate decreases as the absolute value of 
hopping distance increases. This hopping-distance de- 
pendence of the hopping rate was explained quantita- 
tively in two ways: a brick accumulation model and an 
analytical approach. In the brick accumulation model, 
the hopping behavior of the localized Lyapunov vectors 
was explained as the movement of the highest position 
in the brick accumulations. It was shown that using this 
brick model we can calculate the largest Lyapunov ex- 



ponent for quasi-one-dimensional hard-disk systems suc- 
cessfully. On the other hand, in the analytical approach 
the hopping rate was calculated from probability distri- 
butions for brick height differences of nearest neighbor 
sites via multiple summations over possible configura- 
tions that connect two separated highest sites. The result 
is related to the largest Lyapunov exponent. Both of the 
approaches successfully reproduced the hopping-distance 
dependence of the hopping rate for the localized Lya- 
punov vectors of quasi-one-dimensional hard-disk sys- 
tems. 

As a remark, it may be useful to mention a previ- 
ous conjecture for the dynamics of the localized region 
of Lyapunov vectors corresponding to largest Lyapunov 
exponent. Before we started this work, there had been a 
view that the origin of the hopping behavior of localized 
Lyapunov vectors was: 

(Conjecture) The localized region of the Lyapunov vector 
hops to the position of a new grazing collision. 

This conjecture was suggested from the fact that a change 
of Lyapunov vectors in particle collisions may be domi- 
nated by the factor l/(er ■ Ap) in the collision dynamics 
for Lyapunov vector (see Appendix [BJ) , where cr is the 
normalized collision vector and Ap is the momentum dif- 
ference of the colliding particles before the collision. In 
this argument, if two particles collide at a small angle (a 
grazing collision) with a small value of cr ■ Ap, then the 
Lyapunov vector can change significantly and the posi- 
tion of the localized region may move. However, this sce- 
nario cannot be correct for the following reasons. First, 
this conjecture implies that as the position at which a 
grazing collision occurs is random, the hopping distance 
should also be random. This contradicts our numeri- 
cal result that the hopping rate decreases with increas- 
ing hopping distance, as shown in Fig. [21 Second, we 
show in Fig. ^3 the distribution D(cr ■ Ap) of cr ■ Ap 
in general collisions (dashed line), and the distribution 
of collisions causi ng h opping of the localized Lyapunov 
vector (solid line) [48j. It is clear from Fig. EI] that the 
quantity cr ■ Ap is not smaller in collisions which lead to 
jumps of the localized region of Lyapunov vector com- 
pared with the general case. It may also be noted that 
the hopping dynamics of the Lyapunov vector from the 
brick accumulation model described in Sec. IIVI is inde- 
pendent of collision parameters like collision angles, the 
momentum difference of colliding particles, and also the 
quantity cr ■ Ap. This also suggests that the conjecture 
for the origin of the hopping of the localized region of 
the Lyapunov vector cannot be justified. On the other 
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hand, this collision parameter independence of the hop- 
ping behavior in the brick model cannot explain why the 
distribution D(cr ■ Ap) is different from the general case 
and the hopping case in Fig. 1101 This is an open prob- 
lem. 

As another open problem, the fits in Fig. |3| suggest 
that the hopping-distance dependence of the hopping 
rate in the thermodynamic limit seems to be a power 
law: Poo(h) ~ h 13 with (3 1.7. It remains to be deter- 
mined whether this power behavior of the hopping rate 
can be justified analytically in the brick accumulation 
model or the analytical approach discussed in Sec. [V] 

In this paper we showed that the hopping rate of the 
localized region of a Lyapunov vector, and the largest 
Lyapunov exponent, arc well described by the brick ac- 
cumulation model. Then, one may ask a more direct 
numerical check of the justification of the brick accumu- 
lation model in quasi-one-dimensional hard-disk systems, 
for example, to check Eq. J3J or to observe numerically 
an actual brick configuration like that presented in Fig. 
QJ However, such a check of the brick model is not trivial 
for the following reasons. First, the brick accumulation 
model is only justified in the limit of low density, but 
numerical simulations have to be at some finite density. 
This effect appears, for example, as gradual changes of 
the Lyapunov vector component amplitudes, as shown 
in Fig. 1111 which do not appear in the brick accumu- 
lation model. Second, even if we could simulate at an 
extremely low density in which such finite density effects 
can be neglected (although the case presented in this pa- 
per is not at such a low density), the factor — In p in Eq. 

may be too large for an actual numerical calculation. 
Finally, to calculate Lyapunov vectors in this paper we 
used the algorithm developed by Benettin et al [3^, |4jj . 
This algorithm includes intermittent renormalizations of 
Lyapunov vectors, to prevent a divergence of the am- 
plitudes of Lyapunov vectors in numerical calculations, 
but the brick accumulation model does not have such a 
normalization procedure in its dynamics. This difference 
makes a direct numerical check of Eq. J2J difficult in 
hard-disk systems. Different from Eq. itself, the lo- 
calized region of the Lyapunov vectors, which is required 
to calculate the hopping rate, is given simply by parti- 
cle indices with the largest Lyapunov vector component 
amplitude, which is not influenced by such a difference 
of normalization procedure in calculations of Lyapunov 
vectors. 

In this paper, in order to introduce the hopping rate 
we used the property of quasi-one-dimensional systems, 
that the order of particles is an invariant. In this sense, 
it is not trivial to generalize our argument to fully two- 
(or three-) dimensional systems. An effective way to de- 
scribe the dynamics for the localized region of the Lya- 
punov vectors for higher spatial dimensions remains an 
important future problem. Related to this it should still 
be noted that it is known that the clock model version of 
the brick accumulation model itself can be easily general- 
ized to any spatial dimensional case, although in higher 
dimension we do not have the concept of the accumu- 
lation of bricks, as we do for the quasi-one-dimensional 
case. 



Finally, one should notice that the brick accumula- 
tion model (more generally the clock model) used in this 
paper has been justified for hard-disk (or hard-sphere) 
systems only (at least so far). On the other hand, the 
Lyapunov localization is observed not only in many- 
hard-disk systems but also in a wide variety of many- 
body chaotic systems, such as the Kuramoto-Sivashinsky 
model 1221. a random matrix model |24|, map systems 
[23L l25l l44l l45l]. coupled nonlinear oscillators |2fj, etc. 
It should be an important future problem to develop ap- 
proaches to the dynamics of Lyapunov localization in this 
wider class of systems. 
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APPENDIX A: HOPPING RATE OF A 
LOCALIZED LYAPUNOV VECTOR 

In this appendix we give the detailed definition of the 
hopping rate for the localized region of the Lyapunov 
vectors, which is used in this paper. 

As suggested in Ref . [2^, H3 > the normalized Lyapunov 
vectors corresponding to the large Lyapunov exponents 
have non-zero components for only two particles in the 
low density limit, and changes of the non-zero compo- 
nents are caused by particle-particle collisions. Applying 
this characteristic of the Lyapunov vector to the quasi- 
one-dimensional system with periodic boundary condi- 
tions, we can introduce a hopping distance at the 
n-th particle-particle collision as 

h [n] = j n+ i ~ jn ~ N nint j Jn } (Al) 

in the low density limit. Here {j n , j n +l} are the non-zero 
components before the n-th collision and {j n +i, jn+i + 1} 
are the non-zero components just after the n-th collision. 
The function nint jzl is the closest integer to the real 
number x. We assume that changes in the position of 
the localized region of the Lyapunov vector are negli- 
gible during the free-flight interval, so {j n +i > jn+l + 1} 
can also be interpreted as the set of the particle indices 
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FIG. 11: The normalized amplitude 7^ of the Lyapunov vec- 
tor component STj corresponding to the largest Lyapunov 
exponent A' 1 ' as a function of the collision number n t and 
the particle index j in a quasi-one-dimensional system with 
d — 10 4 and N = 50. Here, the collision number interval 
shown in this graph is [k, k + 68] with k — 3600448. On the 
base of this graph is a contour plot of 7) at the level 0.2. 
Three hops of the localized region of the Lyapunov vector are 
visible in this time interval; The first two hops are sharp and 
the last one is gradually. 



whose Lyapunov vector components take non-zero val- 
ues just before the (n + l)-th collision. We count the 
number of times Nx(h) that we see a hop of size h in 
a time-interval T where -[N/2] < h < [N/2]. The 
normalized hopping rate Pjsr{h) can be introduced as 
P N (h)/P N (l) = lim^oo N T (h)/N T (l). 

However, in actual numerical simulations, the particle 
density p is always finite, and non-zero Lyapunov com- 
ponents of more than two particles can often be seen, 
at least down to a density p w 10~ 5 which is the low 
density limit of the numerical simulations in this paper. 
This makes the above definition IjAljl of the hopping dis- 
tance impractical. To explain this point concretely 
we show Fig. 1111 which is a graph of the normalized 
Lyapunov vector component amplitude 7^ defined by 
Eq. as a function of the collision number n t and the 
particle index j in a quasi-one-dimcnsional system with 
N = 50 and d = 10 4 , (namely a density p sa 7.85 x 10~ 5 ). 
In this figure we can recognize three types of hops of the 
localized region of the Lyapunov vector. The first two 
hops keep the non-zero Lyapunov vector components of 
almost two particles sharp enough to apply the definition 
(|A1|) of the hopping distance, but the third hop occurs 
gradually so that no clear hopping time can be deter- 
mined. Note that the localized region of Lyapunov vector 
component amplitudes in the three dimensional plot 1111 
always has a flat top with a width of two-particlcs even 
for the third hop in Fig. I|ll|) . 

In this paper, for concrete calculations, we define the 
localized region of the Lyapunov vector components as 
the particle indices j for which jj > 0.2. (Here, we 
use the similarity between particle positions and parti- 
cle indices given in Fig. ^for the quasi-one-dimensional 
system.) In Fig. ^2 this localized region is approxi- 
mately given by the region surrounded by the contour 
lines (level 0.2) on the base of the graph. We introduce 
the quantity l n as the number of particles satisfying the 

(n) 

inequality 7) > 0.2 just before the n-th collision. Note 

that < 7^ < 1 by definition JTJ of 7^, so l n cannot be 
larger than 5. If l n is always 2 as in the low density limit, 
then we can use the hopping distance definition (|A1|) , but 
in numerical simulations at finite density l n > 2 can hap- 
pen as shown in Fig. 1111 The problem then is how do 
we define the hopping distance /J™' at the n-th collision, 
when l n > 2. 

The definitions of the hopping distance for each case 



FIG. 12: Schematic illustrations of the hopping types for 
the localized region of the Lyapunov vector in quasi-one- 
dimensional systems at low density. The contours represent 
the level 7^ = 0.2. (These are the contours on the base of 
Fig. II 11 1 The vertical dotted line is collision number n at 
which the hopping distance of the localized Lyapunov vec- 
tor is to be determined. The seven illustrations in this figure 
indicate: (a) the case of (l n ,l n +i) — (2,2), (bl) the case of 
(In, In+i) = (2, 3), (b2) the case of (l„, l n+1 ) = (3, 2), (b3) the 
case of (l n ,ln+i) = (3,3), (cl) the case of (l n Jn+i) = (2,4), 
(c2) the case of (l n , ln+i) = (4,2), and (c3) the case of 
(l„,l n +l) = (4,4), where l n is the number of particles in the 
localized region of Lyapunov vector just before the n-th col- 
lision. The numbers in the right-bottom of each illustration 
give the possible values of the hopping distance. 



are categorized as follows. 

(Case (a): l n = 2 — > l n +i = 2 ) Here only two particle 
indices are in the localized region before and af- 
ter the collision, and they are always nearest- 
neighbors, so the hopping distance tv- n ' is given by 
Eq. (23). 

(Case (bl): l n — 2~^ l n +i — 3 ) Here we assume that 
the localized regions of the Lyapunov vector are 
given by {j n ,j n + 1}, {j„+i, j„+i + 1, j n +i + 2} and 
jn+i = in or jn+i = jn ~ 1- We take the value 
of the hopping distance to be 1 (-1) for fv- n ' where 

jn+1 = jn (jn+1 = jn ~ 1)- 

(Case (b2): l n — 3 — > ln+i = 2 ) In this case we assume 
that the localized regions of the Lyapunov vector 
are given by {j n ,jn + 1, jn + 2}, {j n +i,jn+i + 1} 
and j n+ i = j n or j n+x = j n + 1, and the hopping 
distance is in both cases. 

(Case (b3): l n — 3 — * l n +i = 3 ) In this case we assume 
that the localized regions of the Lyapunov vec- 
tor are given by {j n ,jn + 1, j„ + 2}, {j n +ujn+i + 
1) jn+i + 2} and j n+ \ = j n + hi 71 ' with the hop- 
ping distance h^. We take into account the case 
feW = -1,0 or 1 only. 

(Case (cl): l n — 2 — * l n+1 = 4 ) In this case we assume 
that the localized regions of the Lyapunov vec- 
tor are given by {j n ,jn + 1}, {j n +i,3n+i + 1} 
and {jn+i,j'n+i + 1} satisfying j' n+1 = j n and 
br'i+i — jn+i I > 2. The hopping distance is defined 
by Eq. I|A1|1 using these j n +i and j n . 

(Case (c2): l n — 4 — * l n+1 = 2) In this case we assume 
that the localized regions of the Lyapunov vec- 
tor are given by {j„,j„ + 1}, {j„,j« + 1} and 
{jn+ljjn+l + 1} satisfying |j„ - j'J > 2 and 
\jn — jn+i I < 1- The hopping distance hv 1 ' is 
^ n ]=j n+ i-j„ = -l,Qorl. 

(Case (c3): l n = 4 — * l n +i = 4 ) In this case we consider 
only the case in which the localized regions of Lya- 
punov vector are given by {j n ,jn + 1}, {j„,j„ + 1}, 
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FIG. 13: Log-log plots of the normalized hopping rates 
Pjv(/i)/-Pjv(1) for a quasi-one-dimensional system of 50 hard- 
disks as a function of the hopping distance \h\ at different 
densities, d — 10 3 (circles) , 10 4 (triangles) and 10 5 (squares) . 
The density is a function of d given by p = TvR 2 /{(l + d)L 2 y ]. 
The error bars are given by |Pjv(/i) — Pjv(— /i)|/Pjv(1). 



{jn+i,jn+i + 1} and {jn+iJn+l + l } satisfying 

jn+l = jn, j'n+l = j'n and \jn ~ fn\ > 2 - TllC 

hopping distance takes the value in this case. 



(For each case above, the corresponding schematic il- 
lustration is shown in Fig. 1121 ) Here, we use peri- 
odic boundary conditions for the particle index, so that 
the localized regions {N, 1}, {N, 1, 2} and {N - 1, N, 1} 
should be translated to {0,1}, {0,1,2} and {—1,0,1}, 
respectively, in the above definition of the hopping dis- 
tance. In the examples shown in Fig. ^] the first two 
hops of the localized Lyapunov vector can be described 
as case (a), and the third hop is described as cases (cl) 
and (c2). Note that asymmetric definitions of hopping 
distances tv- n ' between cases (cl) and (c2) [and simi- 
larly between cases (bl) and (b2)] are required so that 
we can interpret the hopping distance of the third non- 
zero hop in Fig. 1111 as only —2 in spite of it involving 
both cases (cl) and (c2). Using the above hopping dis- 
tance we count the number Nt{Ii) of hops of distance 
h in time interval T, and introduce the normalized hop- 
ping rate as P N (h)/P N (l) = lim^oo N T (h)/N T (l) for 
h = -[N/2], -[N/2] + 1, • • • , [N/2]. Notice that there 
are possibilities apart from those shown above, but it is 
observed that in numerical simulations the probabilities 
of these are extremely small (for example, more than 96 
percent of all hops could be categorized this way). 

In Fig. 1131 we show the normalized hopping rate 
P/v(/i)/fV(l) in a quasi-one-dimensional system of 50 
hard-disks for d = 10 3 (circles), 10 4 (triangles) and 
10 5 (squares). The hopping rate is almost density- 
independent in this low density range, although it may be 
very slightly larger at the lowest density for \h\ = 2, 3, • • • . 

We also calculate the hopping rate for the brick accu- 
mulation model explained in Sec. II VI in a similar way. In 
the brick accumulation model we can introduce the local- 
ized region of the Lyapunov vectors as the site (particle) 
indices whose brick height are highest. For the brick 
model, cases (bl), (b2), and (b3) above cannot happen, 
and only cases (a), (cl), (c2) and (c3) above are taken 
into account in the numerical calculations. It may be 
noted that in the brick model, case (a) above can hap- 
pen only when the hopping distance is —1, or 1, and 
Eq. (|A1() alone is not enough to calculate the hopping 
distance. This is another reason to take into account 
the case where l n > 2 in the calculation of the hopping 
distance. 



APPENDIX B: CLOCK MODEL FOR 
MANY-HARD-PARTICLE SYSTEMS 



Here we give an extension of the derivation of the clock 
model for the Lyapunov vector dynamics in many-hard- 
particle systems in the limit of low density. We also de- 
rive the formula JSJ for the largest Lyapunov exponent 
from the clock model. The one-dimensional version of the 
clock model is the brick accumulation model used in this 
paper. In particular, we clarify the assumptions needed 
to justify the use of this model to discuss the Lyapunov 
localization. For the basic idea of the clock model, see, 
for example, Refs. [32l l33|. However, note that in this 
appendix we use notation and assumptions that are a 
little different from these references, so that the derived 
clock model is consistent with the discussions in this pa- 
per and able to be compared with the numerical results 
ofRef. 

We consider a D-dimensional system with N hard- 
disks (or hard-spheres, etc.) with identical radius R and 
mass M . We assume that there is no external field in 
the system so that the dynamics is simply free-flights, 
and collisions between two particles. We put 8qj (qj) as 
the spatial part of the Lyapunov vector component (the 
spatial coordinate) of the j-th particle, and Spj (p^) as 
the momentum part of the Lyapunov vector component 
(the momentum) of the j-th particle. 

We take t = t n to be the time of the n-th collision 
which involves particles j n and k n , and define r„ to be 



tn—1 



M 



(Bl) 



so that the n-th free flight time is given by r n M . The 
free flight part of dynamics of the Lyapunov vector is 
represented as 

S qj (t-) = Sq ] (t+_ 1 )+r n Sp J (t+_ 1 ), (B2) 
SPj(t-) = Sp j (t+_ 1 ) (B3) 

Here, the argument refers to the limit of the quantity 
before the n-th collision (— ) or after the n-th collision 
(+). On the other hand, the change in the Lyapunov 
vector in particle-particle collisions is represented as 



Sq jn (t-) + e^5q knjn (t-), (B4) 

Sq kn (t~)-e^Sq kn3n (t-), (B5) 

5 qi (t-), for H{ 3n ,k n }, (B6) 
5p jn (t-) + Q^S Pknjn (t-) 

+Q^Sq knJ Jt-), (B7) 
S Pkn (t-)-Q^Sp knjn (t-) 

-Q^Sq knJ Jt-), (B8) 

Spi(tt) = S Pl (t-), for H{j n ,k n }, (B9) 



Sq t (t+) 

*Pjn(*n) 



where Sq knjn = Sq kn -6q jn , Sp knJn = 5p kn - 8p jn , and 



12 



0["1 and are denned by 



M 



(BIO) 



2i? 



A [n] [n]T- 

x I/- /, * . . ) (Bll) 



with 



ApH = p fcn (t-)-p jB (t-) 



(B12) 
(B13) 



and the (XW) x (VN) identity matrix I. Note that in 
this appendix we introduce all vectors as column vectors, 
so for example, <x["] T Ap[™J is a scalar and cr^Ap^ T is 
a matrix where T is the transpose. For later use we note 
that 

5 Pjn (t+) + SPkM) = Sp jn (t-) + S Pkn (t-) (B14) 



which can be derived from Eqs. (|B7(1 and l|B8(l . For the 
derivation of Eqs. l|B4)l. (|B5| . (|B6| . (|B7| . l|B8|l and JB9|l 

for the Lyapunov vector dynamics, for example, see Ref. 
0. 

We consider a low density case, in which the free flight 
time T n M is large (as the free flight time is inversely 
proportional to the density). This justifies our first ap- 
proximation for the Lyapunov vector: 



<5g, (<J ~ T nSpAt+_ 1 ) 



(B15) 



in the limit of low density, the term containing r„ is much 
larger than the other term on the right-hand side of Eq. 
(|B2ll . This asymptotic relation (|B15() leads to 



S <lk n j n (t n )~ T n 5p k ,{t+_ 1 ) 



(B16) 



from the definition of 5q kn j n (t~) and 8Pk n j n (si-i)- Us- 
ing the relation (jB16|) . Eas." (|B7)l . (JB8J an"d" l|B^) can be 

rewritten as 

+r n Q^S PknJ Jt+_ 1 ), (B17) 
<5p fe j4) ~ S Pk Jt+_ 1 )-e^Sp kn]n (t+_ 1 ) 

-r n Q [n] S PknJn (t+_ 1 ) (B18) 

= <yp,(*»-i). fOT ^{in,M. (Bi9) 



where we have used Eq. I|B3(1 . Note that the spatial part 
of the Lyapunov vector does not appear in the dynam- 
ics described in (|B17|) and HB180 . I|B19|I anymore. The 
first and second terms on the right-hand side of (|B17|) 
and (|B18|) are negligible compared with the third term 
because of the large value of r„ , so we obtain 

6 Pj Jt+) ~ T n Q^6p knj Jt+_ 1 ), (B20) 

SPkM) ~ -TnQ ln] Sp knjn (t+_ 1 ), (B21) 

which lead to 

S Pj Jt+)+Sp kn (t+)~0. (B22) 



On the other hand, for the dynamics given by l|B4|l . I|B5|) 
and (|B6|) for the spatial part of Lyapunov vector we ob- 
tain 



rn [Sp jn (t+_ 1 ) + e^Sp kajn (t+_ 1 )\ (B23) 
rn [SPkJt^) - e^5p knjn (t+_ 1 )] (B24) 
r n S Pl (t+_ 1 ) 1 for lt{j n ,kn}, (B25) 
using Eqs. (|B15|) and (|B16|) . Now we note 

5 Pj n {tn-l) = — o 



+ 



2 

tPjAtn) +^Pk,S t n) , tiPj„k„(tt-l) 



2 2 

_ s Pjn (t+) + spkM) sp^M-i) 

2 2 
~ Jj^lMA (B26) 

where we used Eqs. (|B3|I . (|F314I and (|B22fl . Similarly we 
have 



Spi^jAtt-i) 



(B27) 



Therefore Eqs. I|B23|I and l|B24(l can be rewritten as 

SQjM) ~ ^ [nl ^ fenin (t+_i), (B28) 
<W#) ~ -r n n^S PknJ Jt+_ 1 ) (B29) 
where 17 ^ is defined by 



17 w = -li + e [ 



(B30) 



The a symp totic equations JBT9J>, (|B20|) . (|B21|l. (fB25|) . 
(1B28|) and (|B29|) give the Lyapunov vector dynamics in 
the limit of low density, ft should also be noted that the 
assumptions used to derive this dynamics breaks some 
conservation laws in the original dynamics, for example, 
the quantity $Pj(t) ^ s conserved in the original dy- 

namics (|B3f) . (|B7I) . (|B8|) and [|B9|) but this cannot be 
guaranteed exactly in the low density dynamics ljB19f) . 
(|B20|) and p32l|) . 

Now we consider the Lyapunov vector ST correspond- 
ing to a positive Lyapunov exponent. The positivity of 
the Lyapunov exponent means that the amplitude \5r\ 
of Lyapunov vector diverges exponentially in time. The 
dyna mics given by (|B19|) . (|B2l) . (|B2T|) . (|B25l) . (|B28|) 
and 1B29|I shows that the divergence of \ST\ must come 
from the Lyapunov vector components corresponding to 
colliding particles, as the other Lyapunov vector compo- 
nents diverge at most linearly. For this reason we ne- 
glect the change of the Lyapunov vector components for 
non-colliding particles. Under this assumption the Lya- 
punov vector dynamics for the Lyapunov vector compo- 
nent STj = (SqpSp j) T for the j-th particle is summa- 
rized as 



-6r kn (t+)~ Tn 62l n -V (B31) 
5r l (t+_ 1 ), for H{ Jn ,k n }, (B32) 
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where Sa^ n 1 1 is defined by 
<5H I™ -1 ! = 



n [n] tp knjn (4-i) 



Q [Bl *Pk BiB (d) 



(B33) 



It is essential to note that the Lyapunov vector compo- 
nents 6rj n (tn) and 6rk„(tn) for the colliding particles 
have the same amplitude, and Sa^ n ~^ is independent of 

We assume that the ratio 



A* 



\8pi 



(B34) 



between the amplitudes of the spatial and momentum 
parts of the Lyapunov vector for the i-th particle, are in- 
dependent of the particle index j j^a|. Ref. j^l] suggests 
that the ratio p need not be of order 1 in general. From 
Eq. (|B34|) we have 



1 



V 71 + M' 



(B35) 



Using this assumption we estimate the magnitude of the 
vector (5H[ ,l_1 ] as 



= v /|fi[«]«$ P;Wn (i+-i)| 2 + 

= max{|5p jn (t+_i)| , |^Pfc n (4-i)|} 





f- |g[» 


e [n] 


2 




\$r kn 


(*t 




nN e N 2 - 


f Ql" 


eN 


2 


H 


-M 2 







where is defined by 

[»] = <^Pfc„(*n-l) -<^Pj„(*n-l) 

max { 1 5p kn (t+_x)| , |£Pj n (*n_i)|}' 
and satisfies the inequality 



(B36) 



(B37) 



< 





<5p 


fe„ (*n-l) 


+ 


l*Pi„(*n-l) 




max{ \5p kn (t+ 






*Pi„(Ci) 


} 



< 2. (B38) 



From Eq. (|B36(I we can estimate the magni- 
tudes \Sr jn (t+)\ and through \Sr jn (t+)\ = 

\6r kn (t+)\=T n \53^-% 

Now, we introduce the clock value ICj(n) of the j-th 
particle just after the n-th collision as 



tCj(n) = 
or equivalently, 

\srM)\ 



1, M 

lnp 1^(0)1 ' 



Kj(n) 



1^(0)1 



(B39) 



(B40) 



leading to Eq. 10} . Here p is the particle density whose 
value is < p < 1, and §Tj{Q) is the Lyapunov vec- 
tor component for the j-th. particle at the initial time. 



We choose the initial Lyapunov vector Sr(0) so that its 
components Sr 1 (0),Sr 2 (0),- ■ ■ ,MV_i(0) and Sr N (0) 
have the same order of magnitude, and a larger clock 
value means larger magnitude of the Lyapunov vector 
component, |5-Tj(i n )| . 

We have already used the fact that the free flight time 
increases as the density p decreases. To write down the 
clock value more meaningfully, we use the specific rela- 
tion between the free flight time and the density, namely 
that the free-flight time r n M is approximately inversely 
proportional to the density p at low density; 



(B41) 



Sn/P 



where s n is independent of the density. Using Eq. I|B41|I , 
the expression (|B39|) for the clock value K,j(n) can be 
rewritten as 



KLi(n) 



max{/C in (n- l),K kn {n- 1)} + 1 + A$I 
for I = j n or I = k n 



(B42) 



/C/(n- 1) for I <£{j n ,k n }, 
where A$[™1 is defined by 



1 + M 2 



(B43) 



A$H = — m< St 
mp 



and where we have used l|B31|) . (|B36|) . I|B39|I and 

|(5r , Jre (0)| s» \Sr kn (0)|. Our final assumption to derive 
the clock model is that 



lim A$ l 



0. 



(B44) 



To justify the assumption (|B44|I note that the right-hand 
side of Eq. (|B43I) for the quantity A$[™1 has a factor 
1/lnp which goes to zero in the limit as p — ► 0, and s n , 
and fiW are almost independent of the density p, 
and the magnitude of the vector is finite even in the 
limit p — > from the inequality (|B38|) 46]. Eqs. (IB42I) 
and (|B44|) lead to the clock dynamics 



K,i{n) 



max {JCj n (n — 1), /C kn (n — 1)} + 1 
for I = j n or I = k n 

Ki(n-1) for I £{j n , fc„}, 



(B45) 



in the low density limit, which is closed by the clock 
value Ki{n) only. From Eq. (|B45|) . the dynamics for the 
clock model is expressed as (i) the clock value is changed 
only when the corresponding particle collides, and (ii) the 
clock values of colliding particles are tuned to the same 
value given by 1 plus the larger of the two clock values 
of the particles just before the collision. 

In the quasi-one-dimensional system with periodic 
boundary conditions, particle indices of the colliding par- 
ticles can be taken so that k n — j n + 1 (note that index 
N + 1 is equivalent to 1). Therefore we obtain the dy- 
namics © for the brick accumulation model explained 
in Sec. IIVI as the one-dimensional version of the clock 
model. Moreover, in the one-dimensional case, ICj(n) 
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can be interpreted as the brick height of the j-site just 
after the n-th brick is dropped. 

Finally we derive the formula JSJ for the largest Lya- 
punov exponent from the clock value K-j (n) . In the limit 
of low density p « 1, the amplitude |<$.T(t+)| of the 
Lyapunov vector can be approximated by 



if k: 



\ 



N 

E 



\sr 3 {ti)\ 2 



l(n) 



|^(0)| (B46) 



noting that the sum l^j^n )| 2 i s dominated by the 

Lyapunov vector component amplitude \5rj(t+)\ with 
the largest clock value iO max \n): 

£[mo*]( n ) = max {/Ci(n),/C 2 (n),--- ,£jv(n)} (B47) 

at n, because of the huge factor 1/p. Here, a n is the num- 
ber of particle with the largest clock value K.^ max ^(n), and 
we have also used our assumption that |<5-Tj(0)| is almost 
independent of the particle index j. Now we assume the 



approximate relation 

mol l(n) n - 



1 N 



nN 



(B48) 



3=1 



in the limit of large n. (We have checked this relation 
numerically for the brick accumulation model.) Using 
the relations (IB46I) . (IB48t and t ~ nr we obtain 



A = Urn -ln|jr(t)| 

lim J-^[" laa; ](n)lnp 



lim — !— V/C,(n)lnp (B49) 



3 = 1 



for the Lyapunov exponent A corresponding to the Lya- 
punov vector ST. Therefore we obtain the Eq. JSJ, which 
is independent of the system shape and the number of 
spatial dimensions. 
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we added the fitting parameter a in the fitting function 
Y^,j = i K>j ( n ) = a + fi n in Fig. |S|for consistency, although 
it is almost zero in this case. 
[48] More exactly, the distribution D(er ■ Ap) for non-zero 
hopping cases in Fig. 1101 was calculated as follows. First 
we chose the case satisfying the conditions: (i) the parti- 
cle indices satisfying 7^ > 0.2 are given by {j n ,jn + 1} 
just before the n-th collision and {jn+i,jn+i + 1} just 
after the n-th collision, such that j n 7^ j n +i, and (ii) 
the jn+i-th particle and (j n +i + l)-th particle collide at 
the n-th collision. (Here, if the particle indices for the 
localized region of the Lyapunov vector satisfy 7^ > 0.2 
is {1,7V}, then the particle indices pair {j n ,jn + 1} or 
{j n+ i,j n+ i + 1} above should be taken as {0, 1}, based 
on periodic boundary conditions in the longitudinal di- 
rection of the quasi-one-dimensional system.) Under this 
circumstance we calculate the distribution D(cr ■ Ap) 
from the value of a ■ Ap in the n-th collision between 
these jn+i-th and (j n +i + l)-th particles. 
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